Numerical study of relaxation in electron glasses 
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We perform a numerical simulation of energy relaxation in three-dimensional electron glasses in 
the strongly localized regime at finite temperatures. We consider systems with no interactions, 
with long-range Coulomb interactions and with short-range interactions, obtaining a power law 
relaxation with an exponent of 0.15, which is independent of the parameters of the problem and of 
the type of interaction. At very long times, we always find an exponential regime whose characteristic 
time strongly depends on temperature, system size, interaction type and localization radius. We 
extrapolate the longest relaxation time to macroscopic sizes and, for interacting samples, obtain 
values much larger than the measuring time. We finally study the number of electrons participating 
in the relaxation processes of very low energy configurations. 



I. INTRODUCTION 

Strongly localized systems are characterized by very 
slow relaxation rates due to the exponential dependence 
of the transition rates on hopping length (l]-[3) . For a wide 
range of parameters, the typical times involved are much 
larger than the experimental times and a glassy behav- 
ior is observed. Ben-Chorin et al. [jjj reported on non- 
ergodic transport in Anderson localized films of indium 
oxide and ascribed the phenomena to the hopping trans- 
port in non-equilibrium states. Ovadyahu and Pollak 
performed further experiments on this system that 
clearly demostrate the glassy nature of Anderson insula- 
tors. Glassy behavior may be obtained independently 
of the strength of interactions and regardless of their 
long or short range. In systems with localized states, 
long hopping lengths result in very long relaxation times. 
However, it is thought that there are specific features 
of the glassy relaxation behavior that indeed depend on 
the type and strength of the interactions involved. If 
so, relaxation experiments could be an adequate tool for 
studying the strength of interactions. There has been no 
systematic study of the effects of interactions on the re- 
laxation properties of strongly localized systems, and in 
this paper we try to fill this gap as much as possible. 

Most properties of systems with localized electronic 
states strongly depend on interactions. This is especially 
true for Coulomb glasses where interactions are of a long 
range character. The non-equilibrium properties of these 
systems are affected by dynamic correlations in the mo- 
tion of electrons 0. One-particle densities of states or 
excitations are not enough to encompass the whole prob- 
lem. To deal with such problems, methods were devel- 
oped []|-f7j to obtain the low lying states and energies of 
electron glasses. The states of the system, their energies 
and the transition rates between them constitute the in- 
formation needed to compute non-equilibrium properties. 
We use this information to study energy relaxation for 
systems with no interactions, with long-range Coulomb 
interactions and with short-range interactions. 



In the next section, we describe the model and the 
numerical procedure used. In section III, we study the 
temporal dependence of energy relaxation and, in sec- 
tion IV, we calculate the largest relaxation time T2 and 
its dependence on size and temperature. Finally, in sec- 
tion V, we present results about the number of electrons 
participating in low- energy relaxation processes. 



II. MODEL AND NUMERICAL PROCEDURE 

We consider three-dimensional systems in the strongly 
localized regime, in which quantum overlap energies, h, 
arising from tunnelling are much smaller than the other 
important energies in the problem and are taken into 
account only to the lowest contributing order, i.e., to zero 
order for energies and to first order for transition rates. 
Spin is neglected since exchange energies are proportional 
to t 2 . We use the standard tight -binding Coulomb gap 
Hamiltonian H: 



H = eirii + riinjVij 



(1) 



where ti is the random site energy chosen from a box dis- 
tribution with interval [-W/2, W/2}. For non-interacting 
systems Vij — 0, while Vij = \jr for systems with 
Coulomb interactions and Vij — (0.7/r) 4 is the poten- 
tial chosen for short range interactions. The large value 
of the Hubbard energy is accounted for by disallowing 
double occupation of sites. 

We study systems with sizes from 248 to 900 sites 
placed at random (for short range interactions we only 
consider systems sizes up to 465 sites), but with a min- 
imum separation between them, which we choose to be 
0.5/o where Iq = (47riV/3)^ 1 / 3 and N is the concentration 
of sites. We take e 2 /Iq as our unit of energy and lo as 
our unit of distance. We choose the number of electrons 
to be equal to half the number of sites. We use cyclic 
boundary conditions. 

We use two different numerical algorithms to obtain 
the ground state and the lowest energy many-particle 
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configurations of the systems up to a certain energy. For 
short range interactions, we employ an algorithm that re- 
laxes the system through certain simultaneous n-electron 
transitions ||. The procedure is repeated for different 
initial random configurations of the charges until the con- 
figuration of lowest energy is found ten times. The config- 
urations thus generated were memorized in terms of site 
occupation numbers and of energy, whenever this was less 
than the highest energy configuration in memory storage. 
We complete the set of low-energy configurations by gen- 
erating all the states that differ by one- or two-electron 
transitions from any configuration stored. 

For long range interactions, we use an algorithm that 
consists of finding the low-energy many-particle con- 
figurations by means of a three-step algorithm pOl. 
This comprises local search , thermal cycling fl2f | 
and construction of "neighbouring" states by local re- 
arrangements of the charges [IllJjfl . The efficiency of this 
algorithm is discussed in Ref. [flO|| . In the first step, an ini- 
tial set, S, of metastable low-energy many-particle states 
is created. We start from states chosen at random. These 
states are relaxed by a local search algorithm which en- 
sures stability with respect to excitations from one to 
four sites. In the second step, this set S is improved by 
means of the thermal cycling method, which combines 
the Metropolis and local search algorithms. Lastly, the 
third step completes the set S by systematical investiga- 
tions of the surroundings of the states previously found. 

The transition rate ujjj between configurations / and 
J is taken to be 



Uu = — cxp 1-2 rij/a ) exp I 



kT 



(2) 



for Ej > Ej, and without the second exponential for 
Ej < Ei. In this equation, To is the inverse phonon 
frequency, of the order of 10~ 13 s, a is the localization 
radius, which we take equal to 0.3^o, and X) r u is the 
minimized sum of the hopping lengths of the electrons 
participating in the transition. 

The relaxation process is governed by the master equa- 
tion, which in first order can be written in matrix form 
as p(< + St) = Aip(t), where p is the vector of occupa- 
tion probabilities in the configuration space, and Ai the 
matrix of transition probabilities between states during 
a time, St, given by 



{M)jI ~ { 1 - Ek^i <»iK&t for I = J. 



for I ^ J, 



(3) 



We assume that the system initially occupies a set, 
/C, of to configurations with equal probabilities, that is, 

P'k = V m f° r K e ^-i an d Pl* 1 = for all other L. 
The time evolution of p is governed by the eigenvalues 
Xi and right eigenvectors fa of Ai . We will assume that 
the A, are arranged in decreasing order. Rewriting p( ' 



as a linear combination of the fa, the probability vector 
after n time steps p^™^ is given by 



p (n) = aifa + a 2 fa\^ + a 3 fa\% + 



(4) 



where a, is the i-th component of p(°) in the basis jVi j. 

At long times (large n), Eq. (||) approaches equilibrium 
with time dependences given by A™. Thus, the relaxation 
times are given by 



(5) 



|lnA. t 



in units of St. The final state is = exp(—EM/kT)/Z 
for all M, where Em is the energy of state M, and Z is the 
partition function. Clearly p(°°) is a right eigenvector of 
Ai with eigenvalue 1, since Aip^ 00 ^ — p 1 - 00 * 1 . All the other 
eigenvalues of Ai are smaller than 1, since otherwise the 
system would not tend to the stationary probability dis- 
tribution. The second largest eigenvalue corresponds to 
the largest relaxation time of the system. The addition 
of the other eigenvectors to fa = p(°°', transfers p from 
high energy states to low energy states at various rates. 

We have developed a renormalization method to be 
able to properly handle the huge range of transition rates 
involved. Large values of Tj correspond to A; with values 
which are very close to unity, Eq. ([5]) , and a direct calcu- 
lation of Ti, in units of St, is strongly limited by the nu- 
merical precision of the computer. In order to minimize 
errors, we must choose a St which is as large as possible, 
although this soon yields negative diagonal elements of 
Ai. We overcome this problem using a renormalization 
procedure that allows us to increase St and to simultane- 
ously keep all terms of Ai positive. This procedure forms 
groups of configurations. Each group is made up of con- 
figurations connected between themselves by transition 
rates which are larger than a critical one. The groups 
are clusters in local equilibrium for times greater than 
the inverse of the critical transition rate. Firstly, we take 
a critical transition rate u c . Then for each ujij larger 
than uj c we define a new equilibrium state, M, and sub- 
stitute the original configurations, / and J, by this new 
state, M. The transition rates between M and any other 
configuration K (K ^ I, J) are defined as: 



UKM — WKI + UK J 

UlK , Ujk 

OJMK 



1 + Rm l+R 



AI 



where Rm is given by: 



Rm - — = exp {(Ej - E.,)/k B T} . 

U)JI 



(6) 
(7) 



(8) 



The diagonal matrix elements ujmm are again equal to 1 
minus the sum of the non-diagonal elements of the col- 
umn M multiplied by At. 
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After the matrix Ai has been renormalized by the pre- 
vious procedure, we can increase the time scale to a larger 
interval S't — l/ui c . With this S't we calculate the new 
elements of A4. The eigenvalues of the transition matrix 
will be given now in units of S't(> St). We have checked 
the validity of our renormalization procedure with sev- 
eral samples of small systems where errors are not crit- 
ical. The method minimizes computer errors in the so- 
lution of the eigenproblem as the matrix becomes less 
ill-conditioned, and allows us to consider large systems, 
with matrix elements that differ by many orders of mag- 
nitude. 

III. TEMPORAL DEPENDENCE 

We calculate the temporal dependence of the energy 
of the system when it relaxes from an initial set of high 
energy configurations. At very long times, the longest 
relaxation process involved predominates and we see an 
exponential relaxation. For shorter times, there is an 
almost continuous sequence of relaxation times, which 
gives rise to a power law relaxation (E — E CCi ) oc t~ a . To 
obtain the exponent of this law it is convenient to rep- 
resent the absolute derivative of the energy with respect 
to time. In Fig. [l] we show |di?/dt| versus time (in units 
of To) in a double logio plot for a sample with Coulomb 
interactions and 248 sites. The continuous curve cor- 
responds to a temperature T = 0.004, and the dashed 
curve to T = 0.005. The straight line is a fit to the data 
in the non-exponential part of both curves, and its slope 
is equal to —1.15. So the power-law exponent for relax- 
ation is a = 0.15. This exponent is basically independent 
of temperature for all the systems considered. 
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FIG. 1. Double logio plot of the temporal derivative of 
the relaxation energy versus time for a system with Coulomb 
interaction, for T = 0.004 (solid curve) and 0.005 (dashed 
curve). The straight line corresponds to power-law relax- 
ation, and has a slope equal to —1.15. t is given in units of 

TO- 



We have also studied energy relaxation for systems 
with short-range interactions and for non-interacting 
systems. The results for short-range interactions are very 
similar to those for Coulomb interactions. The power- 
law exponent is roughly 0.15 and the largest relaxation 
time is of the same order of magnitude as for Coulomb 
systems. In Fig. || we show |di?/d£| as a function of time 
in a double logio plot for a non-interacting system with 
N = 248 sites. The continuous curve is for T = 0.004, 
and the dashed curve for T = 0.005. The slope of the 
straight line is again equal to —1.15. There are two dif- 
ferences between the results for interacting and for non- 
interacting systems. The longest relaxation times are 
shorter for the latter, and the power-law regime is not 
very well defined in the absence of interactions. Both 
figures give the rate of relaxation |d-E/di| at any time. 
At very small t, the interacting systems relax faster than 
the non-interacting systems. A possible explanation of 
this is that in the excited state of the interacting systems 
some electrons get very close to each other. In the initial 
stages of relaxation these electrons hop away from elec- 
trons in the nearest neighbors sites, the whole process 
being very fast. 
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FIG. 2. Double logio plot of |<LE/dt| versus time for a 
non-interacting system, for T — 0.004 (solid curve) and 0.005 
(dashed curve). The straight line has a slope equal to —1.15. 
t is given in units of to. 

Several samples have been checked and in all of them 
we obtain similar results to Figs. |] and ||. Two features 
characterize our relaxation process, the exponent a of the 
power-law regime and the longest relaxation time. The 
exponents a do not appreciably vary from sample to sam- 
ple, nor with temperature or with the type of interaction. 
On the other hand, the longest relaxation time drastically 
changes from sample to sample and with changes in tem- 
perature, size and the range of interaction. On average, 
this time increases with the size of the system and with 
the strength of the interactions. In the next section we 
study the longest relaxation time in detail. Now we shall 
analyze exponent a. 
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Temporal relaxation can be described as a sum of par- 
allel exponential relaxation processes, each with its own 
different relaxation time, Tj. The energy, E, of the sys- 
tem can be written as a function of time, t, as follows 

E(t) = C * eX P I" ) + ^Eq (9) 

i>l V Ti > 

where Cj is the product of the i-th component of the ini- 
tial occupation vector, a,, and the energy associated to 
the eigenvector fa. This energy is the sum of the compo- 
nents of 4>i multiplied by the corresponding energies. Eftq 
is the equilibrium energy, i. e. E-^q = E(t — > oo). In 
Fig. H we plot (ci/n) exp(— i/ri) for the 30 largest eigen- 
values of Ai, excluding Ai = 1, as a function of time for 
a sample with Coulomb interactions and of size N = 465. 
The solid line represents the temporal derivative of the 
actual energy as a function of time. This curve is below, 
but very close to, the envelope of the curves correspond- 
ing to the individual relaxation processes. Note how the 
combination of several simple exponential relaxation pro- 
cesses gives rise to a power law relaxation. 
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FIG. 3. Exponential relaxation processes coming from sev- 
eral eigenvalues of in a double logio plot. A power law 
arises from the combination of all of them. 

Surprisingly, a is fairly independent of temperature, 
size, type of interaction considered, and localization ra- 
dius, facts for which we do not have any interpretation. 
Anyway, the robustness of the exponent could be a sig- 
nature of self-organized criticality. Similar trends have 
been found in experimental measurements of the excess 
conductance of 2D samples excited far from equilibrium 
M. In the absence of magnetic field, the power law ex- 
ponent of these measurements ranges between 0.27 and 
0.29, diminishing with the strength of the magnetic field. 

Our results point to the difficulty in extracting infor- 
mation about the effects of interactions from the power 
law exponent. Nevertheless, the type of interaction sig- 
nificantly affects the longest relaxation times. 



IV. LONGEST RELAXATION TIME 

We also study the longest relaxation time, T2, as a 
function of temperature and the size of the sample for 
systems with Coulomb interactions, with short-range in- 
teractions and for non- interacting systems. In Fig. 4 
we plot (log 10 T2) versus the inverse of the temperature 
for the three types of interactions mentioned, Coulomb 
(solid lines), short-range (dotted-dashed lines) and no- 
interactions (dashed lines). The number of sites consid- 
ered are N = 248, 341, 465, 744 and 899, for long range 
interactions and for non-interacting systems; for short 
range interactions we did not use the two largest sizes. 
T2 increases with sample size, and thus the smallest sam- 
ple corresponds to the lowest curve, and so on. () denotes 
averages over site configurations. Fluctuations in r 2 from 
sample to sample are very large and, as is the case with 
most properties of disordered systems, one has to aver- 
age the common logarithm of T2, rather than t-i itself. 
The curves extend over the range of validity of the re- 
sults. The 'high' temperature limit T max depends on the 
energy range AE spanned by the configurations stored. 
We choose T max = 0.1 AS. The low temperature limit 
arises from the discrete nature of the spectrum of config- 
urations and we take it as being equal to the mean energy 
spacing of the ten lowest energy configurations Ae. 

From Fig. [| we can conclude that the longest relax- 
ation time depends strongly on the type of interaction. 
T2 is one order of magnitude larger for interacting than 
for non-interacting systems. As we will see, this effect is 
much larger when extrapolated to macroscopic sizes. 
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FIG. 4. Average of log 10 T2 versus 1/T for systems with 
Coulomb interaction (solid lines), a short range interac- 
tion (dotted-dashed lines) and systems without interactions 
(dashed lines). The size of the systems stems from 248 sites 
(lowest line) to 900 sites (highest lines). For short range in- 
teractions the largest size considered is 465 sites. 

In order to extrapolate the previous results to macro- 
scopic sizes we plotted (log 10 T 2) as a function of L^ 13 at 
a fixed temperature for different values of the exponent 



4 



f3. L = iV 1 / 3 is the length of the side of the system, and 
N is the number of sites. We found that the results for 
the three types of interactions fit straight lines fairly well 
when /3 = 1. In Fig. || we show log 10 t 2 versus L^ 1 for sys- 
tems with Coulomb interactions (dots), short-range inter- 
actions (diamonds) and without interactions (squares). 
The horizontal dashed line represents a macroscopic time, 
say, one day (« 10 18 to). The temperature chosen in 
this plot is T = 0.0025, which is valid for the four sizes 
employed in both types of interactions. The size of the 
symbols used roughly corresponds to the standard devi- 
ation of log 10 r 2 . The crossing point of each straight line 
with the vertical axis is the extrapolation of t 2 to macro- 



scopic sizes. The results are t^ ^ w 10 31±1 to = 10 1 
s (Coulomb interactions) rs 10 n±1 s (short-range 

interactions) and r^ 1 « 10 5±1 s (no interactions). It 
is clear from this figure that the longest relaxation time 
drastically increases with the strength of interactions, al- 
though these results have to be taken with care as they 
are extracted from a very long extrapolation. 
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sequence of configurations with decreasing energies. For 
each transition at T = 0, the shorter the hopping length, 
the faster the corresponding transition rate. From each 
configuration, the system chooses the nearest one (in 
terms of ^ r) from those with less energy. With this 
in mind, we have computed for all low-energy configura- 
tions the closest one of smaller energy, and have stored 
the number of electrons n participating in the transition. 

In Fig. |^ we show the number of electrons, n, of the 
fastest transition from an initial configuration as a func- 
tion of the number of this configuration for a Coulomb 
interacting sample with 900 sites. At very low energies, 
the relative importance of many-electron transitions in- 
creases. The proportion of transitions with a fixed num- 
ber of electrons greater than one (n > 1) increases with 
decreasing energy. Obviously, in the non-interacting case 
all processes are one-electron transitions. 
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FIG. 6. Number of electrons participating in the fastest 
transition as a function of the order of the initial configura- 
tion. 



FIG. 5. Average of log 10 T2 versus L _1 at T" 1 = 400 for 
non-interacting systems (squares) and systems with Coulomb 
(dots) and short range (diamonds) interactions. 

The results presented in Figs. || and || correspond to a 
localization radius a = 0.3l . For larger values of a, the 
relaxation times will decrease, as can be deduced from 
Eq. (||). We have found empirically that a change in 
a causes a change in r 2 of approximately A log 10 r 2 ~ 
3A(a _1 ). The values of r^ 00 ' are so large for interacting 
systems that we would expect non-ergodic behaviour for 
these systems even for much larger localization radii than 
the one considered here. 

V. VARIABLE NUMBER RELAXATION 

At zero temperature, the relaxation process is down- 
ward in energy and we can assume that the fastest pro- 
cess always dominates, corresponding to a well defined 



VI. CONCLUSIONS 

Our numerical results of relaxation in localized elec- 
tronic systems show a power law behavior with an ex- 
ponent close to 0.15 and independent of all the param- 
eters and type of interactions considered. At very long 
times, we obtain exponential relaxation with a charac- 
teristic time that strongly varies with size, localization 
radius and type of interaction. The extrapolation of this 
characteristic time to macroscopic sizes predicts values 
much larger than the typical experimental times, espe- 
cially for the interacting cases. The strength of interac- 
tions in experiments performed on these systems can be 
deduced from their longer relaxation times. 
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